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ABSTRACT 

Sharp density features in protoplanetary discs, for instance at the edge of a magneti- 
cally dead zone, have recently been proposed as effective barriers to slow down or even 
stop the problematically fast migration of planetary cores into their central star. Den- 
sity features on a radial scale approaching the disc vertical scale height might not exist, 
however, since they could be Rayleigh (or more generally Solberg"H0iland) unstable. 
Stability must be checked explicitly in one-dimensional viscous accretion disc models 
because these instabilities are artificially eliminated in the process of reducing the full 
set of axisymmetric equations. The disc thermodynamics, via the entropy stratifica- 
tion, and its vertical structure also influence stability when sharp density features are 
present. We propose the concept of Rayleigh adjustment for viscous disc models: any 
density feature that violates Rayleigh stability (or its generalization) should be dif- 
fused radially by hydrodynamical turbulence on a dynamical time-scale, approaching 
marginal stability in a quasi-continuous manner. 

Key words: Accretion, accretion discs - hydrodynamics ~ instabilities - planetary 
systems: formation - planetary systems: protoplanetary discs. 



1 INTRODUCTION 



In current theories of planet formation, one of th e most accepted and studied scenarios is the core accretion model (see reviews 
by Lissauei 19931 : Paoaloizou Sz Terauem 2006I . and references therein). In this model, a Jovian planet is built by accreting 
gas onto a sol id planetary core. However, before the core can grow t o a critical mass of ~5-15 for runaway gas accretion 
to occur (e.g. IPoUack et al]ll996l : lAlibert et aPbood : iRafikovlEooel ). type I migration due to angula r momentum exchange 
with the natal protoplanetary gas disc can drive the core into the host star i n a short time-scale (e.g. iGoldreich fc Tremaine 



198' 



Ward 



1997 



Tanaka. Takeuchi fc Wardlboo j iMenou fc Goodman 



(e.g. iHartmann et al 



19981: 



Sicilia-Aguilar et al. 



2001 



i 



Hillenbranc 



di boosi ) 



2004 ) compared to a typical disc lifetime of <10 Myr 



A proposed mechanism to withstand rapid mward type 1 migration is to induce a steep radial density gradient m the ga s 
disc (jMatsumura fc Pudritzlboosl . bood : iMasset et allbooel : iMatsumura. Pudritz fc Thommesll2007l . l20od : llda fc Linlboosi ). 
In particular, since the gas inside the orbit of a planetary core exerts a positive torque on the core whereas the gas outside 
exerts a negative torque, a large negative radial densit y gradient in the gas can generate a sufficiently positive torque to halt 
type I migration. Matsumura fc Pudritz 1 20051 . 20061 ) propose that the transition region in a protoplanetary disc between 
hydromagnetic turbulence and dead zone can locally provide such a density barrier. On the other hand, by including the 

' ^) suggest that a positive torque on 



Masset et 



al 



(2001 



co-rotation torque exerted by the gas in the vicinity of a protoplanet 
the protoplanet can be generated by a positive local den sity gradient. Ida fc LinI ( 20081 
density bump around the ice line to act as a barrier, and lSchlaufman. Lin fc Idal ( 2009 
mock orbital configurations of extraso lar planetary systems. 

However. iPapaloizou fc LinI (|1984| ) have already pointed out that any local variation occurring on the order of a disc scale 
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height may render the disc Rayleigh unstable. The sharp feature may then be washed out by the resulting hydrodynamical 
turbulence. Given that a large local density gradient is necessary to impact type I migration, it is thus critical to examine the 
stability of such a feature. 

We examine the Rayleigh stability of local density features in ^ In ^ we introduce the concept of Rayleigh adjustment 
as a possible solution to deal with such features in one-dimensional viscous disc models. In §4] we discuss the need to generalize 
these results by accounting for the possibly important role of entropy stratification and the vertical structure. We conclude 
with some general comments in ^ 



2 RAYLEIGH STABILITY 



Since magnetically-coupled regions in a protoplanetary disc are unstable to the magneto-rotational instability (jBalbus fc Hawlev 
19981) and unlik ely to develop sharp features, our main concern is the hydrodynamical stability of a magnetically dead zone 
i Gammielliggil l. The Rayleigh stability criterion is a necessary and su fficient condition for the local axisymmetric stability 
of an inviscid differentially rotating fiuid system l|Chandrasekhar^ll961^ . It states that the specific angular momentum must 
monotonically increase with cylindrical distance R from the central rotation axis in a fiow for strict stability: 

^(^^0) > 0, (1) 

where v^, is the azimuthal velocity of the fiow. Near the mid-plane of a protoplanetary gas disc, which is supported by the 
central stellar gravity and the local pressure gradient, force balance in the radial direction requires that 
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V4, 



2 RdP 
p oR 



(2) 



where vk is the Keplerian velocity at R, p is the gas density, P is the gas pressure, and any radial infall has been neglected. 
Substituting equation ([2]) into equation ([T| gives the stability criterion 
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(3) 



where we have used vk oc R~^^^. 

If variations in pressure, i.e. density or temperature for an ideal gas law, occur on a length scale of order i? in a disc, the 
second term in equation @ is of the order c^, where Cs — y^'yP/p is the sound speed and 7 is the adiabatic index (typically 
7/5 for diatomic H2 in the protoplanetary context). In a protoplanetary disc, it is usually the case that vk ^ Cs and thus the 
disc is Rayleigh stable irrespective of the sign of the variation. However, if any variation occurs on a length scale of order the 
disc scale height, H = R(cs/vk), the same term becomes of the order (caR/H)^ ~ vj^, assuming vertical hydrostatic balance 



holds as usual (e.g. Frank. King fc Raine 20021 '). In other words, the second term in equation ([Sj is of the same order as the 
first term. Therefore, the disc is prone to Rayleigh instability and sensitive to the exact profile of pressure variations. 

Without loss of generality (but see SI) ^6 adopt an ideal gas law for the equation of state and focus our stability analysis 
on the disc mid-plane. We also assume for simplicity that the sound speed Cs is slowly-varying with cylindrical distance R, as 
compared to the density p: \{dcs/dR)/{dp/dR)\ <S Cs/pE|Then, to leading order, the stability condition (eq. [3]) approximately 
becomes 

7?2 d^p P2 '^"^2 

n2 



V Cs 



+ 



R^ (dp\\md^^ 

p 9i?2 p2 \q^) + p dR~ 

For density variations on a length scale ~ H, the fourth term in equation Q is of the order R/H • 



(4) 

Vk I Cs and thus may 
Therefore, both 



be neglected. On the other hand, both the second and the third terms are of the order (R/H)^ ~ {vr/cs 
the steepness and the curvature of the density profile are important to determine the stability of the disc. 

Since the third term in equation Q is always negative, both steeply ascending and descending density profiles may make 
the disc Rayleigh unstable. For example, we consider a two-fold density drop near R — Ro of the form 

[3 1 , , fR- Ro^ 



p{R) = po 



1 

- tanh 
2 



^ R-Ro \^ 



(5) 



where po is the density after the drop while e controls the width of the drop. At R — Rg, the first derivative of the density 
profile reaches maximum while the second derivative vanishes. Substituting equation ((SJ in the stability condition @ and 
evaluating it at = Rq gives (vk/cs)'^ — (i?o/e)^/9 > 0, where we have ignored the fourth term in equation @. Since 
H/R = Cs/vk, it requires that e > H/3 and the maximum slope allowed in this case is \dp/dR\ ~ 3po/2H. This illustrates 
how Rayleigh stability sets an upper limit on the steepness of a density profile. 



^ Our analysis could be generalized to address a steep radial temperature profile but we note that one would expect such a feature to 
be reduced by horizontal radiative transport (see also JJlfor the role of entropy stratification). 
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Furthermore, the second term in equation Q indicates that a disc may be Rayleigh unstable where the second derivative 
of the density profile is negative. In regions where the profile is concave down, the term sets a lower limit to the radius of 
curvature. In particular, a cusp pointing upward cannot be stable due to an infinite negative second derivative. Smoothness 
is a necessary condition in these regions. Locations where pressure or density maxima occur in a protoplanetary disc are of 
particular interest in that solid grains are captured, but Rayleigh instability may also come into question at these locations 
in a nontrivial manner. To illustrate this point, we consider a density bump with a Gaussian form, 



PiR) = Po 



, H 

i H cxp 

e 



(6) 



Superposed on a disc of uniform mid-plane density po, the bump encloses a constant mass (on the order of the mass in an 
annulus of mid-plane density po and width H at R — Rq) and has a width proportional to e, provided that e <g Rq. The 
most negative second derivative occurs at the peak, where the first derivative vanishes. Evaluating the stability condition (0} 
at R^ Ro results in (vk/cs) - 2H{Ro/ey /{H + e) > 0. Since H/R = c^/vk, the minimum width allowed for stability is 
e ^ H. Therefore, a local density maximum cannot be Rayleigh stable if the width of the bump is less than about a disc scale 
height. 



3 RAYLEIGH ADJUSTMENT 



Once a certain region in a protoplanetary disc reaches a state of Rayleigh instability, we conjecture that the ensuing hydrody- 
namical turbulence will redistribute material in the disc in such a way that any narrow feature will be smoothed and broadened 
till the disc recovers a state that is (marginally) Rayleigh stable. The time-scale ta for a disc to recover a Rayleigh-stable state 
might be close to the dynamical time-scale of the disc, i.e. tji > , where f^K is the Keplerian angular velocity, although 
a slower growth rate may be expected as one app roaches marginal stability. By contrast, standard a-model of thin accretion 
discs has a viscous time-scale ti, ^ R^/aHcg (e.g. Frank et al.l 2003 ) . Given that a < 1 and H <^ R, > {R/H)^tR 3> tR. In 



other words, any model of viscously-evolving protoplanetary discs should maintain Rayleigh stability quasi-steadily. 

If the condition tv S> tn holds, it is useful to implement an adjust ment scheme to guarantee that a disc model r emain s 
Rayleigh quasi-stable at all time. In evolutionary models like those of Matsumura et al. 1 2007 . 20091 ) or Ida fc Liiil ( 20081 ) . 
for instance, this kind of adjustment could be applied at any time-step whenever Rayleigh instability is detected, before one 
proceeds with the migration torque calculation. The purpose of the adjustment is to capture the final state of the adjustment 
process when the disc just restores its Rayleigh stability and rotational equilibrium (implied by eq. [2]), without any detailed 
knowledge of the Rayleigh hydrodynamical turbulence occurring on small temporal or spatial scales. 

In the hydrodynamical context of interest here, Rayleigh adjustment must involve some level of mass and angular mo- 
mentum redistribution in the disc. To the extent that the disc reaches a new state of rotational equilibrium (described by 
eq. [2]) after Rayleigh adjustment, however, the adjustment itself can be described as a transition from one state of rotational 
equilibrium to another such state, under the action of mass redistribution. Therefore, this process may be phenomenologically 
described by the radial mass diffusion equation 



dt RdR\ dRj' 



(7) 



where E is the vertically-integrated column density and D is the mass diffusion coefficient. We use E instead of p under 
the hypothesis that the transport is vertically global. Whether or not the disc evolution due to Rayleigh turbulence can 
be described by fast radial mass diffusion remains to be demonstrated by full hydrodynamical simulations. In this regard, 
equation ([7| can only be deemed to be an approximation to the actual adjustment process. 

We hereby describe a simple algorithm to implement Rayleigh adjustment, using the diffusion equation ([7]). As a first 
step, we ignore any spatial variations and assume that f is a constant to focus our discussion on how to include Rayleigh 
adjustment into viscously-evolving disc models. The equation can be discretized as 



Ef 



^ DAt 



AR 



(e;+i-2e," + e7_,) + ^(e 



(8) 



where E" is the density at ii = Rj and step n, At is the time increment of the adjustment scheme, and AR — Rj+i — Rj 
is the grid spacing. We have adopted a regular mesh in equation (8) for simplicity. From von Neumann stability analysis. 
At < AR^/2D. Substituting At = AR^/2D in equation gives 



= ^ E 



■j+i 



+ 47?, 



(9) 



A Rayleigh unstable density profile can then be treated as an initial guess and equation ([9} can be iterated to the point when 
the density first becomes everywhere Rayleigh stable, under the constraint that rotational equilibrium is satisfied throughout 
(eq. [2]). This corresponds to the completion of the turbulent redistribution process induced by Rayleigh instability. Note that 
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Figure 1. Initially Rayleigh unstable {solid linos) and final Rayleigh-stable {dashed and dot-dashed lines) profiles for (a) a density drop 
(eq. [5]) and (b) a density bump (eq. [6]) processed through the Rayleigh adjustment scheme discussed in With the constant diffusion 
coefficient adopted in this algorithm, the final relaxed density profiles do not depend on the width of the initially unstable profiles. 



the magnitude of the diffusion coefficient D plays no role in determining the final state of the Rayleigh adjustment, which can 
be seen in equation (|9} , since we have chosen the maximum time-step possible to minimize the number of iterations needed to 
find the adjusted state. We have verified that the final profile does not depend on the time-step, as long as the von Neumann 
stability condition is satisfied, and that the total mass of the disc is conserved. 

We emphasize that since the diffusion time-scale of this process is about L^/D, where L is the characteristic length 
scale of the density profile, sharp Rayleigh-unstable density features (with L <^ R) will diffuse much more rapidly than any 
density structure on large scales. Therefore, in first approximation, the global disc structure will remain largely unaffected 
by this Rayleigh adjustment even for a globally constant D. Note that such a quasi-instantaneous adjustment scheme can 
only be justified if the process leading to Rayleigh instability (e.g. viscous mass pile-up) operates much more slowly than the 
adjustment itself. 

As a demonstration, we have implemented this adjustment scheme on the two model density profiles, equations ([Sjl 
and ([6|, previously discussed in 321 For concreteness, we adopt Ro = 1 and H — 0.1. The parameter e is chosen such that the 
density profile is Rayleigh unstable. The profile is discretized with 1000 grid points in the range 0.5 ^ R ^ 2 and diffused via 
equation ([9|. As shown in Fig. [l] profiles of different initial widths relax to the same enlarged width at marginal Rayleigh 
stability, which is roughly one disc scale height, H. Even though these results are based on a constant mass diffusion coefficient, 
they are consistent with the expectation that the width of any density profile is intrinsically limited to ~ _ff by Rayleigh 
stability. One expects any variation narrower than '-^ _ff to be smoothed and broadened, although details of this process 
may depend on the specific Rayleigh adjustment prescription adopted. We note that, in addition to redistributing mass and 
angular momentum, localized hydrody namical turb ulence in the disc could also directly influence the migration process (e.g. 



angular momentum, locaiizea nyaroay namicai turp uience m tne aisc couia also girectiy mnuence tne migration process {e.g. 
Laughhn. Steinacker fc Adamsll2004l : iNelson 20055; I Johnson. Goodman fc Menoul [2OO6I : lOishi. Mac Low fc MenoulboO?! ). an 



effect which is not captured by our adjustment scheme. Potential consequences for the disc thermodynamics have also been 
ignored. 

We note that the mass diffusion equation ^ employed here conserves mass in the disc but it does not conserve angular 
momentum during the adjustment process, under the assumption that rotational equilibrium holds before and after Rayleigh 
adjustment, as specified by equation ([2)|. At first, this would seem to be a serious limitation of the simple mass redistribution 
scheme proposed here. One may be tempted to improve such a scheme by using instead a more standard viscous equation for 
mass and angular momentum transport in a thin disc, of the form 



1/2 A/„v??i/2n 



i^--(.Ei^ 



(10) 



dt ^ RdR 

where u is the effective disc viscosity (e.g. Frank et al. 2002[ ). However, even though this equation does conserve mass and 



angular momentum in a Keplerian disc, it would be as inadequate as equation (O if used in a Rayleigh adjustment scheme. 
Indeed, equation (|10p . or its generalizations for non-Keplerian rotation laws, assumes that the disc angular velocity profile is 
fixed in time. It is precisely the breakdown of this assumption during Rayleigh adjustment, dv^/dt 7^ 0, going from one state 
of rotational equilibrium to another such equilibrium, which prevents this entire class of equations from accurately tracking 
angular momentum conservation during a Rayleigh adjustment episode. To the extent that one is interested in describing the 
global evolution of discs in rotational equilibrium with approximate equations such as equation (|10|l . one may then choose 
to consider the small violations in angular momentum conservation occurring during Rayleigh adjustments as an acceptable 
source of errors, given that it is an unspecified detail at the level of approximation of the viscous theory. This, in our view. 
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motivates the use of the simplest possible mass diflFusion scheme, under the constraint of rotational equilibrium, as we have 
proposed here with equation ((7)|. 

To summarize, the general outline of a Rayleigh adjustment procedure would be as follows. At each time-step of a one- 
dimensional evolutionary disc model, evaluate the stability of the disc profile using equation ([3} or a generalization (see 311 
below). If the profile is stable, move on to the next time-step of the evolution. Otherwise, use the unstable profile as the 
initial condition and iterate equation ((O)! , or a more sophisticated form of equation ((TJ , till the profile has relaxed to marginal 
stability. Then use the relaxed profile to resume the evolutionary model. Two important assumptions underlying this scheme 
are (i) that the agent driving the disc viscous evolution acts on a time-scale much longer than the disc dynamical time, (ii) that 
the relaxed density profile is independent of the pre-adjusted profile, as illustrated in Fig. [1] for the simple scheme proposed 
here. 



4 ROLE OF ENTROPY STRATIFICATION AND VERTICAL STRUCTURE 

So far, our discussion has focused exclusively on the Rayleigh stability of a purely radial stratification of angular momentum 
in the disc mid-plane. More generally, entropy stratification will also be present in the protoplanetary disc and contribute to 
its axisymmetric stability if sharp density features exist. The full generalization of Rayleigh's stability analysis in the presence 
of entropy st ratification leads to the two necessary and sufficient Solberg-H0iland criteria for axisymmetric stability (e.g. 
Tassoullll97i) : 



^i + ^i + ^^jT^>o, (11) 



dR 

dP \ f 1 dR'^n^ din Pp-'' 1 dR'^n^ din Pp- 
"dzj [ R^ OR dZ R^ dZ dR~ 



where 



1 dP din Pp-'' 1 dP din Pp~ 



-ypdR OR ' ^ -ypdZ dZ 



> 0, (12) 



(13) 



Q = Vtj,/ R IS the angular velocity, and Z is the vertical coordinate along the rotation axis. The term A'^^ -I- is the sum of 
the cylindrical radial and vertical components of the squared Brunt-Vaisala frequency. 

Equations (llip - (|12[) make it clear that the stratifications of entropy and angular momentum both contribute to axisym- 
metric stability. Since these criteria were derived with respect to any linear axisymmetric perturbation in the {R, Z) plane, 
they address the cylindrical radial and the vertical components of both stratifications. In the presence of a sharp radial density 
feature, the radial component of the entropy stratification, which is normally negligible in a disc, can make a significant con- 
tribution to local stability. To illustrate this possibility, let us consider a case again where the sound speed is slowly- varying 
with cylindrical distance R as compared to the density, and let us restrict perturbations to be purely cylindrical radial so that 
vertical stratification plays no role (djdZ ^ in eqs. [ll)-[13p. In this limit, axisymmetric stability reduces to the much 
simpler criterion 

' ^^'"%0, (14) 



" i?3 dR 

which is the equivalent of equation ([TJ with an additional contribution from radial entropy stratification. According to the 
same scaling analysis as used in iJUfor equation Q, in the presence of a radial density gradient with length scale i/, neglecting 
any radial temperature gradient on such scales. 



r2 ^. <^B ( . 1 \ ^. "If 1 



The magnitude of the radial entropy stratification, as measured by A'^^, could thus be comparable to that of the angular 
momentum stratification and contribute to the radial stability of the disc configuration. 

This suggests that, in the presence of sharp density gradients, it may no longer be possible to ignore the disc thermo- 
dynamics for stability considerations. Perhaps even more importantly, equations (|lip - (|12|l show that the disc radial stability 
can no longer be studied independently of its vertical structure. In the presence of radial density features with length scale 
~ H , radial and vertical gradients become comparable in magnitude. The disc stability is then determined by a detailed 
balance between the radial and vertical stratifications of angular momentum and entropy (eqs. [Il]-[l2]). This makes the 
stability analysis significantly more complicated than considered so far, especially in the case of dead zones with a poorly 
known vertical structure. It may also imply that a more elaborate treatment than the simple Rayleigh adjustment scheme 
proposed in 33] is required for unstable discs. 
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5 DISCUSSION AND CONCLUSION 



The possibility that dis cs with sharp radial density features may be subject to Rayleigh instability has already been noted 
by various authors (e.g. Papaloizou fc Lin 1984 : Matsumura et al. 2007 . 20091 ) . We have reconsidered this issue here and have 
suggested that stability be checked explicitly in viscous evolutionary disc models which develop such features. In viscous disc 
models, the full set of axisymmetric equations is reduced to a single radial diffusion equation of the type shown in equation (|10p . 
Typically, the equation reduction is achieved by imposing a steady angular velocity profile in the disc, n(7?, t) — fl{R) (often 
chosen to be Keplerian, like in eq. [10]), and by neglecting the radial acceleration term in the radial momentum equation. 
These are simplifications which contribute to eliminating the class of axisymmetric instabilities discussed here. If sharp density 
features are to be considered as viable solutions t o slow down or even stop protoplanetary migration, with many import ant 
consequences for planet formation scenarios (e.g. iMatsumura et all booj I2OO9I : llda fc Lin|[2008l : ISchlaufman et al.ll2009l '). a 
careful examination of stability in viscously evolved models seems warranted. 

The va rious arguments we hav e put forward are not meant to imply that the s harp density features considered in various 
works (e.g. Matsumura et al. 2007 . 20091 : Ida fc Lin 20081 : Schlaufman et al. 20091 ) are axisymmetrically unstable. However, 
they suggest that an approximate scaling analysis based on a representative length-scale ~ H may not be sufficiently accurate 
to evaluate the axisymmetric stability of a disc, which depends on the detailed shape of the density profile. In addition, 
since both the vertical structure and the entropy stratification of the disc can influence its stability when sharp features are 
present, a careful treatment of the disc thermodynamics may be required. The extent to which density features can or cannot 
grow to be as sharp as a disc scale height ma y significantly affect planetary migr ation because differential Lindblad torques 



H (e.g. 



Ward 1997: Matsumura et al 



2007) and co-rotatio n torques are very sensitive 



are also applied over a length-scale ^ ^ ^ 

to the local density i Masset et al. 20061 ) and thermodynamic ( Paardekooper fc Papaloizou 20081 ) conditions. For those discs 



which develop unstable density profiles as a result of slow viscous evolution, for instance mass pile-up in a magnetically dead 
zone, we have proposed that Rayleigh adjustment schemes be implemented. Although more detailed work would be needed to 
incorporate such a scheme within the framework of a viscous disc solver, our goal is to suggest that this type of adjustment 
schemes can be used to evolve more consistently discs with sharp features on viscous evolutionary time-scales. 

Finally, we remark that our discussion has been restricted to the stability of axisymmetric perturbations in protoplanetary 
discs. It has been suggested that dis cs with narrow density features are also susceptible to non-axisymmetric instabilities 



b een su gg 

( Li et al.ll200ol . l200ll : lLvra et al.ll2009h . To the extent that such instabilities lead to mass diffusion on a fast dynamical time- 



scale, it should be possible to model their effects on disc evolutionary time-scales with an adjustment scheme similar to the 
one proposed here. 
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